{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Supply chain experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipyleaflet as ipyl\n",
    "import geopandas as gpd\n",
    "import shapely.wkb \n",
    "import requests\n",
    "import urllib\n",
    "import os\n",
    "import ee\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utils\n",
    "**df_from_carto**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_from_carto(account, query):\n",
    "    \"\"\"\n",
    "    It gets data by querying a carto table and converts it into a GeoDataFrame.\n",
    "    \"\"\"\n",
    "    urlCarto = f\"https://{account}.carto.com/api/v2/sql\"\n",
    "    \n",
    "    sql = {\"q\": query}\n",
    "    r = requests.get(urlCarto, params=sql)\n",
    "    \n",
    "    data = r.json()\n",
    "    \n",
    "    df = gpd.GeoDataFrame(data.get(\"rows\"))\n",
    "    if 'the_geom' in df.columns:\n",
    "        # Change geometry from WKB to WKT format\n",
    "        df['geometry'] = df.apply(lambda x: shapely.wkb.loads(x['the_geom'],hex=True), axis=1 )\n",
    "        df.drop(columns='the_geom', inplace=True)\n",
    "        if 'the_geom_webmercator' in df.columns:\n",
    "            df.drop(columns=['the_geom_webmercator'], inplace=True)\n",
    "        df.crs = {'init': 'epsg:4326'}\n",
    "        df = df.to_crs({'init': 'epsg:4326'})\n",
    "        \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basemaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mongabay_Paper = {\n",
    "    \"url\": 'https://api.mapbox.com/styles/v1/mongabay/ckai1bze50hmj1ipohkupl0qo/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoibW9uZ2FiYXkiLCJhIjoiY2s4N25oZG82MGFyejNsb2lnN3ZrbG1jbyJ9.coW3do99wi_PmnW6OylFbA',\n",
    "    \"name\": 'Mongabay_Paper'\n",
    "}\n",
    "\n",
    "topography_url = 'https://api.mapbox.com/styles/v1/casius/ck801p48x1hsd1iqwe67w1lac/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiY2FzaXVzIiwiYSI6ImJDMkpucTQifQ.5rm4_TsT8_PH8TzOY2V3FQ'\n",
    "topography_basemap = ipyl.TileLayer(\n",
    "    url=topography_url,\n",
    "    layers='topography',\n",
    "    format='image/png',\n",
    "    name='Topography',\n",
    "    opacity=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cerrado boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ikersanchez/anaconda3/envs/geoenv/lib/python3.8/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n",
      "  return _prepare_from_string(\" \".join(pjargs))\n"
     ]
    }
   ],
   "source": [
    "account = 'p2cs-sei'\n",
    "query =('SELECT * FROM brazil_biomes WHERE name = \\'CERRADO\\'')\n",
    "df = df_from_carto(account, query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_file('../data/cerrado/cerrado.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "locations = []\n",
    "geometry = df['geometry'].iloc[0]\n",
    "for polygon in geometry:  \n",
    "    x, y = polygon.exterior.coords.xy\n",
    "    locations.append(list(zip(list(y), list(x))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sentinel 2 Cloud Free Composite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CloudMaskS2(image):\n",
    "    \"\"\"\n",
    "    European Space Agency (ESA) clouds from 'QA60', i.e. Quality Assessment band at 60m\n",
    "    parsed by Nick Clinton\n",
    "    \"\"\"\n",
    "    AerosolsBands = ['B1']\n",
    "    VIBands = ['B2', 'B3', 'B4']\n",
    "    RedBands = ['B5', 'B6', 'B7', 'B8A']\n",
    "    NIRBands = ['B8']\n",
    "    SWIRBands = ['B11', 'B12']\n",
    "\n",
    "    qa = image.select('QA60')\n",
    "\n",
    "    # Bits 10 and 11 are clouds and cirrus, respectively.\n",
    "    cloudBitMask = int(2**10)\n",
    "    cirrusBitMask = int(2**11)\n",
    "\n",
    "    # Both flags set to zero indicates clear conditions.\n",
    "    mask = qa.bitwiseAnd(cloudBitMask).eq(0).And(\\\n",
    "            qa.bitwiseAnd(cirrusBitMask).eq(0))\n",
    "\n",
    "    return image.updateMask(mask)\n",
    "\n",
    "def CloudFreeCompositeS2(startDate, stopDate):\n",
    "    ## Define your collection\n",
    "    collection = ee.ImageCollection('COPERNICUS/S2')\n",
    "\n",
    "    ## Filter \n",
    "    collection = collection.filterDate(startDate,stopDate)\\\n",
    "            .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))\\\n",
    "            .map(CloudMaskS2)\n",
    "\n",
    "    ## Composite\n",
    "    composite = collection.median()\n",
    "    \n",
    "    ## normDiff bands\n",
    "    normDiff_band_names = ['ndvi', 'ndwi']\n",
    "    for nB, normDiff_band in enumerate([['B8','B4'], ['B8','B3']]):\n",
    "        image_nd = composite.normalizedDifference(normDiff_band).rename(normDiff_band_names[nB])\n",
    "        composite = ee.Image.cat([composite, image_nd])\n",
    "    \n",
    "    return composite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Landsat 8 Cloud Free Composite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "def maskL8sr(image):  \n",
    "    # Bits 3 and 5 are cloud shadow and cloud, respectively.\n",
    "    # Get the pixel QA band.\n",
    "    qa = image.select('pixel_qa');\n",
    "    # Both flags should be set to zero, indicating clear conditions.\n",
    "    mask = qa.bitwiseAnd(1 << 3).eq(0).And(qa.bitwiseAnd(1 << 5).eq(0))\n",
    "    return image.updateMask(mask)\n",
    "\n",
    "def CloudFreeCompositeL8(startDate, stopDate):\n",
    "    ## Define your collection\n",
    "    collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')\n",
    "\n",
    "    ## Filter \n",
    "    collection = collection.filterDate(startDate,stopDate).map(maskL8sr)\n",
    "\n",
    "    ## Composite\n",
    "    composite = collection.median()\n",
    "    \n",
    "    return composite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SRTM Digital Elevation Data 30m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = ee.Image('USGS/SRTMGL1_003')\n",
    "elevation_30m = dataset.select('elevation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## USGS National Elevation Dataset 10m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = ee.Image('USGS/NED');\n",
    "elevation_10m = dataset.select('elevation');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Area of Interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon = {\n",
    "  \"type\": \"FeatureCollection\",\n",
    "  \"features\": [\n",
    "    {\n",
    "      \"type\": \"Feature\",\n",
    "      \"properties\": {},\n",
    "      \"geometry\": {\n",
    "        \"type\": \"Polygon\",\n",
    "        \"coordinates\": [\n",
    "          [\n",
    "            [\n",
    "              -46.79351806640625,\n",
    "              -12.149430892248033\n",
    "            ],\n",
    "            [\n",
    "              -45.274658203125,\n",
    "              -12.149430892248033\n",
    "            ],\n",
    "            [\n",
    "              -45.274658203125,\n",
    "              -11.372338792141125\n",
    "            ],\n",
    "            [\n",
    "              -46.79351806640625,\n",
    "              -11.372338792141125\n",
    "            ],\n",
    "            [\n",
    "              -46.79351806640625,\n",
    "              -12.149430892248033\n",
    "            ]\n",
    "          ]\n",
    "        ]\n",
    "      }\n",
    "    }\n",
    "  ]\n",
    "}\n",
    "\n",
    "polygon = {\n",
    "  \"type\": \"FeatureCollection\",\n",
    "  \"features\": [\n",
    "    {\n",
    "      \"type\": \"Feature\",\n",
    "      \"properties\": {},\n",
    "      \"geometry\": {\n",
    "        \"type\": \"Polygon\",\n",
    "        \"coordinates\": [\n",
    "          [\n",
    "            [\n",
    "              -105.35682678222656,\n",
    "              39.98501230042296\n",
    "            ],\n",
    "            [\n",
    "              -105.1944351196289,\n",
    "              39.98501230042296\n",
    "            ],\n",
    "            [\n",
    "              -105.1944351196289,\n",
    "              40.056789440203374\n",
    "            ],\n",
    "            [\n",
    "              -105.35682678222656,\n",
    "              40.056789440203374\n",
    "            ],\n",
    "            [\n",
    "              -105.35682678222656,\n",
    "              39.98501230042296\n",
    "            ]\n",
    "          ]\n",
    "        ]\n",
    "      }\n",
    "    }\n",
    "  ]\n",
    "}\n",
    "\n",
    "polygon = {\n",
    "  \"type\": \"FeatureCollection\",\n",
    "  \"features\": [\n",
    "    {\n",
    "      \"type\": \"Feature\",\n",
    "      \"properties\": {},\n",
    "      \"geometry\": {\n",
    "        \"type\": \"Polygon\",\n",
    "        \"coordinates\": [\n",
    "          [\n",
    "            [\n",
    "              -105.34412384033203,\n",
    "              39.95554342883535\n",
    "            ],\n",
    "            [\n",
    "              -105.2493667602539,\n",
    "              39.95554342883535\n",
    "            ],\n",
    "            [\n",
    "              -105.2493667602539,\n",
    "              40.065723429743045\n",
    "            ],\n",
    "            [\n",
    "              -105.34412384033203,\n",
    "              40.065723429743045\n",
    "            ],\n",
    "            [\n",
    "              -105.34412384033203,\n",
    "              39.95554342883535\n",
    "            ]\n",
    "          ]\n",
    "        ]\n",
    "      }\n",
    "    }\n",
    "  ]\n",
    "}\n",
    "\n",
    "region = polygon.get('features')[0].get('geometry').get('coordinates')\n",
    "roi = ee.Geometry.Polygon(region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Composites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2_composite = CloudFreeCompositeS2(startDate='2019-01-01', stopDate='2019-12-31')\n",
    "s2_composite = s2_composite.clip(roi)\n",
    "\n",
    "viz = {'bands': ['B8', 'B4', 'B3'], 'max':3000, 'gamma': 1}\n",
    "s2_NIR = s2_composite.visualize(**viz)\n",
    "\n",
    "viz = {'bands': ['B4', 'B3', 'B2'], 'max':3000, 'gamma': 1}\n",
    "s2_RGB = s2_composite.visualize(**viz)\n",
    "\n",
    "viz = {'min':-0.75,'max':0.75, 'bands':['ndvi']}#, 'palette': ['#306466','#9cab68', '#cccc66', '#9c8448', '#6e462c']}\n",
    "s2_NDVI = composite.visualize(**viz)\n",
    "\n",
    "viz = {'min':-0.75,'max':0.75, 'bands':['ndwi']}#, 'palette': ['#8bc4f9','#c9995c', '#c7d270', '#8add60', '#097210']}\n",
    "s2_NDWI = composite.visualize(**viz)\n",
    "\n",
    "L8_composite = CloudFreeCompositeL8(startDate='2019-01-01', stopDate='2019-12-31')\n",
    "L8_composite = L8_composite.clip(roi)\n",
    "\n",
    "viz = {'min':0,'max':3000, 'gamma':1.4, 'bands':['B4', 'B3', 'B2']}\n",
    "l8_RGB = L8_composite.visualize(**viz)\n",
    "\n",
    "viz =  {'min': 1250, 'max': 2250, 'gamma': 1.6} #{'min': 0, 'max': 1000}\n",
    "elevation_30m_image = elevation_30m.clip(roi).visualize(**viz)\n",
    "\n",
    "viz =  {'min': 1250, 'max': 2250, 'gamma': 1.6} #{'min': 0, 'max': 1000}\n",
    "elevation_10m_image = elevation_10m.clip(roi).visualize(**viz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1a26b6a0087540c4b8e77ba3e11e7ef7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[40.010628254521095, -105.29674530029338], controls=(ZoomControl(options=['position', 'zoom_in_text…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "images = {'Elevation 30 m': elevation_30m_image,\n",
    "          'Elevation 10 m': elevation_10m_image,\n",
    "          'S2_RGB': s2_RGB,\n",
    "          'L8_RGB': l8_RGB}\n",
    "\n",
    "ee_tiles = 'https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}' \n",
    "\n",
    "m = ipyl.Map(layers=(ipyl.basemap_to_tiles(Mongabay_Paper),),center = (roi.centroid().getInfo().get('coordinates')[1], roi.centroid().getInfo().get('coordinates')[0]), zoom=12)\n",
    "\n",
    "m.add_layer(topography_basemap)\n",
    "\n",
    "for name, image in images.items():\n",
    "    mapid = image.getMapId()\n",
    "    tiles_url = ee_tiles.format(**mapid)\n",
    "    \n",
    "    tile_image = ipyl.TileLayer(\n",
    "        url=tiles_url,\n",
    "        format='image/png',\n",
    "        name=name,\n",
    "        opacity=1\n",
    "    )\n",
    "    \n",
    "    m.add_layer(tile_image)\n",
    "\n",
    "for n, location in enumerate(locations):\n",
    "    polygon = ipyl.Polygon(locations=location, color=\"black\", weight=2, fill=False, name='Countries')\n",
    "    m.add_layer(polygon)\n",
    "\n",
    "control = ipyl.LayersControl(position='topright')\n",
    "m.add_control(control)\n",
    "m.add_control(ipyl.FullScreenControl())\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Export images to drive**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = {'Elevation_10m': elevation_10m_image,\n",
    "          'S2_RGB': s2_RGB,}\n",
    "\n",
    "for name, image in images.items():\n",
    "    task = ee.batch.Export.image.toDrive(\n",
    "        image = image,\n",
    "        description = f'export_{name}_to_Drive', \n",
    "        folder='Mongabay',\n",
    "        fileNamePrefix = name,\n",
    "        scale = 10, \n",
    "        region = roi.getInfo().get('coordinates'),\n",
    "        maxPixels = 1e13,\n",
    "        fileFormat = 'GeoTIFF'\n",
    "    )\n",
    "    task.start()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Export images as `png`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "visSave = {'dimensions': 1024, 'format': 'png', 'crs': 'EPSG:3857', 'region':region}\n",
    "path = '../data/'\n",
    "\n",
    "images = {'Elevation_30m_2': elevation_30m_image,\n",
    "          'Elevation_10m_2': elevation_10m_image,\n",
    "          'S2_RGB_2': s2_RGB,\n",
    "          'L8_RGB_2': l8_RGB}\n",
    "\n",
    "for name, image in images.items():\n",
    "    url = image.getThumbURL(visSave)\n",
    "    \n",
    "    if os.path.exists(path+f'{name}.png'):\n",
    "        os.remove(path+f'{name}.png')\n",
    "    \n",
    "    urllib.request.urlretrieve(url, path+f'{name}.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
